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ABSTRACT 

Direct imaging observations have revealed spiral structures in protoplanetary disks. Previous studies 
have suggested that planet-induced spiral arms cannot explain some of these spiral patterns, due to 
the large pitch angle and high contrast of the spiral arms in observations. We have carried out three 
dimensional (3-D) hydro dynamical simulations to study spiral wakes/shocks excited by young planets. 
We find that, in contrast with linear theory, the pitch angle of spiral arms does depend on the planet 
mass, which can be explained by the non-linear density wave theory. A secondary (or even a tertiary) 
spiral arm, especially for inner arms, is also excited by a massive planet. With a more massive planet 
in the disk, the excited spiral arms have larger pitch angle and the separation between the primary and 
secondary arms in the azimuthal direction is also larger. We also find that although the arms in the 
outer disk do not exhibit much vertical motion, the inner arms have significant vertical motion, which 
boosts the density perturbation at the disk atmosphere. Combining hydro dynamical models with 
Monte-Carlo radiative transfer calculations, we find that the inner spiral arms are considerably more 
prominent in synthetic near-IR images using full 3-D hydrodynamical models than images based on 
2-D models assuming vertical hydrostatic equilibrium, indicating the need to model observations with 
full 3-D hydrodynamics. Overall, companion-induced spiral arms not only pinpoint the companion’s 
position but also provide three independent ways (pitch angle, separation between two arms, and 
contrast of arms) to constrain the companion’s mass. 

Subject headings: accretion, accretion disks - planet-disk interaction - protoplanetary disks - stars: 
protostars 


1. INTRODUCTION 

Recent high-resolution direct imaging observations 
have revealed spiral structure in three protoplanetary 
disks around Herbig Ae/Be stars: SAO 206462 (Muto 
et al. 2012; Garufi et al. 2013), MWC 758 (Grady et 
al. 2013; Benisty et al. 2015), and HD 100546 (Currie et 
al. 2014). The polarized intensity has been measured in 
these observations to gain higher contrast between the 
disk and the central star. While the thermal emission 
from the central star is unpolarized, the scattered light 
from the disk is polarized. In these near-infrared (near- 
IR) polarized intensity images, two spiral arms with 
roughly 180^ rotational symmetry are present in both 
SAO 206462 and MWC 758, similar to the grand design 
in a spiral galaxy (e.g. the Whirlpool Galaxy M51). The 
spiral arms also exhibit a high contrast against the back¬ 
ground disk. The polarized intensity of the spiral arm 
is several times higher than that of the region outside 
the spiral arm. It should also be noted that, since the 
dust scattering opacity is quite large, these observations 
only probe structure high up at the disk atmosphere (e.g. 
several disk scale heights) where the last dust scattering 
surface is. 

In addition to spiral patterns, these three disks also 
have gaps or holes (SAO 206462 with a submillimeter 
cavity of 46 AU, MWG 758 with a cavity of 73 AU from 
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Andrews et al. 2011, and HD 100546 with a cavity of 10 
AU revealed by SED fitting from Bouwman et al. 2003), 
which indicates that they are members of the protoplan¬ 
etary disk class called transitional disks (Espaillat et al. 

2014) . One scenario to explain both spiral patterns and 
gaps/holes is that these disks harbor low-mass compan¬ 
ions (e.g. young planets) which can open gaps and excite 
spiral waves at the same time (e.g., Baruteau et al. 2014). 

However, there are two difficulties in explaining the ob¬ 
served spiral patterns using planet-induced spiral wakes. 
Eirst, the large pitch angle of all the observed spiral arms 
suggests that the disk has a relative high temperature 
(e.g. ^250 K at 70 AU for MWC 758, Benisty et al. 

2015) . In linear theory, spiral waves are basically sound 
waves in disks, and the pitch angle of the spiral arms 
is directly related to the sound speed in the disk. Us¬ 
ing the linear theory, the best fit models for both SAO 
206462 (Muto et al. 2012) and MWC 758 (Grady et al. 
2013; Benisty et al. 2015) suggest that the disk aspect 
ratio {H/R with = Cg/U) at R ^ 100 AU is around 
0.2 which is too large for any realistic disk structure. Eor 
example, even if the stellar irradiation is perpendicular 
to the disk surface R the maximum disk temperature due 
to the stellar irradiation is cfT{R)^ = L^j^irR^ so that 

H/R (X oc . Assuming a 2 Mq central star with 
10 1/0 luminosity, the maximum temperature is ^70 K 

at 100 AU and i7/R is only ^0.1. Since H/R (x it 
is very difficult to make H/R ^0.2. 

Second, the observed spiral arms exhibit much higher 
brightness contrasts than suggested by the synthetic ob- 

^ In reality, the stellar irradiation to the disk is not that efficient 
since the light from the star impinges very obliquely on the disk. 
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servations based on two dimensional (2-D) planet-disk 
simulations. Juhasz et al. (2015) have calculated the po¬ 
larized scattered light images by combining 2-D hydro- 
dynamical simulations with 3-D Monte-Carlo radiative 
transfer (MCRT) simulations. Vertical hydrostatic equi¬ 
librium has been assumed to extend the 2-D simulation 
to the third dimension (the vertical direction). They find 
that a relative change of about 3.5 on the spiral arms in 
the surface density is required for the spirals to be de¬ 
tectable. This value is a factor of eight higher than what 
is seen in their hydro dynamical simulations. 

In this paper, we first point out that the pitch angle 
formula derived from the linear theory, which has been 
used in almost all previous spiral arm modeling efforts, 
does not apply to the high planet mass cases. Spiral 
wakes that are excited by high mass planets (e.g. 1 Mj) 
become spiral shocks which propagate at speeds faster 
than the local sound speed (Goodman & Rafikov 2001, 
Rafikov 2002 ). Thus, the pitch angle difficulty above 
can be alleviated by considering the non-linear exten¬ 
sion of the spiral shock theory. We also show that spiral 
arms (especially the inner arms) have complicated non¬ 
hydrostatic 3-D structure. Such structure can lead to 
strong density perturbation at the disk surface resulting 
in a corrugated shape of its atmosphere. Since near-IR 
observations only probe the shape of disk surface, this 
effect alleviates the second difficulty mentioned above. 
In Dong et al. (2015), we have combined MCRT simu¬ 
lations with hydrodynamical simulations from this pa¬ 
per and demonstrated that planet-induced inner spiral 
arms can explain recent near-IR direct imaging observa¬ 
tions of SAG 206462 and MWC 758. We note that, since 
the planets that we have proposed are outside the spi¬ 
ral arms, we cannot explain the gaps discovered at small 
radii in these transitional disks. Other mechanisms, e.g. 
another planet or photoevaporation, are needed to ex¬ 
plain these gaps. 

Before we introduce our numerical method in §3, we 
provide the theoretical background in §2. The shape of 
the spiral wakes will be studied in §4, and their 3-D struc¬ 
ture will be presented in §5. After a short discussion in 
§ 6 , we summarize our results in §7. 

2. THEORETICAL BACKGROUND 

As a result of planet-disk interaction, a spiral arm 
forms due to the constructive interference of density 
wakes with different azimuthal wavenumbers m excited 
by the planet at Lindblad resonances. In the linear den¬ 
sity wave theory, the m-th Fourier component of the 
planet potential excites the density wave having m spiral 
arms 

5{R, (/), t) = k^(R)dR+m{4,-n^t)] 

where 5 is any perturbed quantity associated with the 
wave, 5q{R) is its complex amplitude, kR{R) is the ra¬ 
dial wave vector, and ftp is the planet orbital frequency. 
Thus, the wave has the same phase along the curve sat¬ 
isfying dR/dcj) = —m/kR^R). The pitch angle (/3) of 
the equal phase curve satisfies tan/3 = \dR/{Rd(j))\^ so 
/d=tan“^|m/[/ci^(R)R] |. Using the dispersion relation¬ 
ship for density waves in the large m limit and far from 
the launching point, rin?{fl{R) — ftpY ~ clkR{R)‘^^ we 
have ^=tdin~^[cs/{R\ ft(R) — ftp\)]. Because P is indepen¬ 
dent of m, different m modes can constructively interfere 


to form the one armed spiral wake (Ogilvie & Lubow 
2002). If the equal phase curve is integrated from the 
planet’s position (R^, 0 ^), the shape of the wake far from 
Rp is given by Rafikov ( 2002 ) and Muto et al. ( 2012 ) as 

tip 

_ 

Rp J ^ f T ^ f 

I I 

1 + 77 1 — a -h rj 

where hp = R/R is the disk aspect ratio at Rp, ft(R) oc 
R~'^, and the sound speed Cs{R) oc R~^. 

However, when the planet is massive enough, the above 
linear density wave theory breaks down. Linear waves 
excited by planets will steepen to shocks (Goodman & 
Rafikov 2001, Rafikov 2002 , Dong, Rafikov & Stone 2012 , 
Duffell & MacFadyen 2012 , Zhu et al. 2013) after they 
propagate over a distance 
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where 7 is the adiabatic index, and Mth is the disk ther¬ 
mal mass 


^th 




( 4 ) 


When Mp > Mth^ the spiral waves will immediately be¬ 
come spiral shocks after they are excited around the 
planet. Unlike the linear wake which follows Equation 
the spiral shock will expand away from the trajec¬ 
tory predicted by Equation Q. Thus, if there is a mas¬ 
sive planet in the disk, using Equation ([^ to fit the shape 
of the spiral shocks will predict an incorrect disk aspect 
ratio and temperature. 


3. NUMERICAL SIMULATIONS 
3.1. Method 

To study density wakes/shocks excited by planets, we 
have carried out both 2-D and 3-D hydrodynamical sim¬ 
ulations using Athena++. Athena++ is a newly de¬ 
veloped grid based code using a higher-order Godunov 
scheme for MHD and the constrained transport (GT) to 
conserve the divergence-free property for magnetic fields 
(Stone et al. , in preparation). But in this paper, we do 
not include magnetic fields and only solve hydrodynam¬ 
ical equations using Athena++. Gompared with its pre¬ 
decessor Athena (Gardiner & Stone 2005, 2008; Stone et 
al. 2008), Athena- 1 — 1 - is highly optimized and uses flexi¬ 
ble grid structures, allowing global numerical simulations 
spanning a large radial range. Eurthermore, the geo¬ 
metric source terms in curvilinear coordinates (e.g. in 
cylindrical and spherical-polar coordinates) are carefully 
implemented so that angular momentum is conserved ex¬ 
actly (to machine precision), which makes the code ideal 
for global disk simulations. 

Our simulations use the adiabatic equation of state 
(EoS) with the adiabatic index 7 = 1 . 4 . A simple orbital 
cooling scheme has been applied to mimic the radiative 
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cooling process in disks. In 3-D simulations, we have 
adopted 

dE _ E CypEj^rprp 

dt tcool 

where p and E are the density and the internal energy per 
unit volume, while in 2-D simulations, we have adopted 


dE 

dt 


E — CyETir 

tcool 


(6) 


where S and E are the disk surface density and the in¬ 
ternal energy per unit area. Cy = k/{pmu{'y — 1)) is the 
heat capacity per unit mass, k is the Boltzmann constant, 
p is the mean molecular weight, and rriu is the atomic 
mass unit. The cooling time tyooi can be written in the 
dimensionless form as Tyooi = tcooi^iR)- We fix Tyooi 
to be a constant in each simulation. With this scheme, 
the disk temperature is relaxed to the background disk 
temperature (Tiyy) determined by stellar irradiation. In 
our simulations, Tiyy is set to be the initial disk temper¬ 
ature. To estimate Tyooi in a realistic disk, we use the 
grey atmosphere approximation (Hubeny 1990) for the 
radiative cooling. 


dE 

dt 



1+t2 ’ 


( 7 ) 


where a is the Stefan-Boltzmann constant, r = {E/2)kr 
is the optical depth in the vertical direction, hZR is the 
Rosseland mean opacity, and Ty is the midplane temper¬ 
ature. Assuming E = CyETy and using Equations § 
and Q, we can derive 

. _ __ 1 + ^^ /ox 

cool lQ^(^T^+Tl^)^T, + Tirr) T ■ 

Approximating the polynomial of Ty and Tiyy in the de¬ 
nominator of Equation with [max(Tc, and as¬ 

suming the central star is 1 Mq, we have 


Tcooi - 0.002 

1+t2 

X- 

T 


/lOOAU 


(60K)3 

[m.ax{Tc,TirrW 

( 9 ) 


Thus, Tcool can vary dramatically at different radii in 
disks. Using the minimum mass solar nebular model, 
Tcool is ^10^ at 1 AU and 10“^ at 100 AU. We have car¬ 
ried out three sets of simulations with Tyyoi = 10“^, 1, 
and 100. They are respectively labeled as ISO, Tl, T2 at 
the end of their names in Table 1. Simulations with fast 
cooling {Tyooi = 10“^) are equivalent to locally isother¬ 
mal simulations. Simulations with Tyooi = 100 are basi¬ 
cally adiabatic simulations considering the timescale of 
the simulations is several tens of orbits, and simulations 
with Tc = 1 should be between isothermal and adiabatic 
simulations. Somewhat surprisingly, we find that simu¬ 
lations with Tyooi = 1 are qualitatively similar to those 
with Tyooi = 100. This similarity suggests that the spiral 
waves in disks with Tyooi = 1 behave adiabatically. We 
think that this is due to the short timescale for the flow 
to move across the spiral wake. The spiral wake has a 
typical width smaller than the disk scale height, and the 
background flow moves across the wake at nearly Kep¬ 
ler ian speed. Thus, when the fluid travels in and out of 


the wave/shock, its response time is much smaller than 
the orbital time and should behave adiabatically even in 
disks with Tyooi = 1- Considering this similarity, in most 
part of the paper, we only show results with Tyooi = 10“^ 
and 1. 


TABLE 1 
Models 


2-D 

Run 

Mp 

Tool 

Domain 

Resolution 


Mj 



R 


Rx (j) 

CMIISO 

0.01 

10-^ 


0.2,10' 


1280x2048 

CMITI 

0.01 

1 


0.2,10^ 


1280x2048 

CMITIOO 

0.01 

100 


0.2,10^ 


1280x2048 

CM2ISO 

1 

10-^ 


0.2,lO' 


1280x2048 

CM2T1 

1 

1 


0.2,10^ 


1280x2048 

CM2T100 

1 

100 


0.2,10^ 


1280x2048 

CM3ISO 

6 

10-^ 


0.2,10' 


1280x2048 

CM3T1 

6 

1 


0.2,10^ 


1280x2048 

CM3T100 

6 

100 


0.2,10^ 


1280x2048 

3-D 

Run 

Mp 

Tool 

Domain 

Resolution 


Mj 



r X 6 


r X 6 X (f) 

STHIN 

0.000316 

10-^ 

[0.5,2]X 

[f-0.2,f+0.2] 

456x128x2048 

SMIISO 

0.01 

10-^ 

[0.3,3]X 

[f-0.6 

,f-h0.6] 

256x128x688 

SMITI 

0.01 

1 

[0.3,3]X 

[f-0.6 

,f-h0.6] 

256x128x688 

SMITIOO 

0.01 

100 

[0.3,3]X 

[f-0.6 

,f-h0.6] 

256x128x688 

SM2ISO 

1 

10-^ 

[0.3,3]X 

[f-0.6 

,f-h0.6] 

256x128x688 

SM2T1 

1 

1 

[0.3,3]X 

[f-0.6 

,|-h0.6] 

256x128x688 

SM2T100 

1 

100 

[0.3,3]X 

[f-0.6 

, 1+0-6] 

256x128x688 

SM3ISO 

6 

10-^ 

[0.3,3]X 

[f-0.6 

,f-H0.6] 

256x128x688 

SM3T1 

6 

1 

[0.3,3] X 

f-0.6 

,f-1-0.6 

256x128x688 

SM3T100 

6 

100 

[0.3,3] X 

[f-0.6, f-1-0.6] 

256x128x688 

SM2IS01ong 

1 

10-'> 

[0.3,3] X 

[f-0.6, f-hO.6] 

256x128x688 

SM2ISOcyl 

1 

10“® 

[0.3,3] 

X [-0.8,0.8] ^ 

296x176x690 ^ 

SM2ISOconstT 

1 

10-5 

[0.3,3] X 

[f-0.6, f-hO.6] 

256x128x688 


^ This is the domain size in the R x Z direction since cylindrical 
coordinates have been used. 

^ The resolution is for R x Z x (j) with the cylindrical coordinate 
system. 


We have also varied the planet mass to be 0.01, 1, and 
6 Mj in our main set of simulations, which are labeled 
as Ml, M2, and M3 in their names respectively. Here, 
we have defined Mj as 0.001 of the central star’s mass. 
The thermal mass (Equation® for the h = 0.1 disk is 
^ Mj. Thus, waves excited oy a 0.01 Mj planet are 
in the linear regime and waves from a 6 Mj planet are 
in the highly non-linear regime. To compare with Eig- 
ure 2 in Tanaka et al. (2002), we have also carried out 
a thin disk simulation with Hp/Rp = 10“^’^ (STHIN 
in Table 1). The thermal mass for such a thin disk is 
only Mfh =0.0316 Mj. Thus, in order to ensure that the 
waves are in the linear regime, we choose the planet mass 
of O.OlMfh = 3.16 X 10“^Mj in this thin disk simulation. 
To avoid the divergence of planet potential, a smoothing 
length of 0.1 Rp has been applied for M2 and M3 cases 
(0.1 Rp is close to the disk scale height and the Hill ra¬ 
dius of the planet). Eor the thin disk case which has a 
very small mass planet, we choose a smoothing length of 
6xlO~^Rp, roughly the length of two grid cells. Eor the 
low mass planet cases (Ml), a smoothing length of 0.02 
Rp, which is also roughly the length of two grid cells in 
these simulations, has been adopted. Planets are fixed 
in circular orbits at R = 1, and the indirect potential, 
which is due to the center of the coordinate system is 
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at the star instead of the center of the mass, has been 
included. We have run the simulations for 10 planetary 
orbits. We choose this timescale because it is longer than 
the sound crossing time throughout the whole disk so 
that density waves/shocks have established, while it is 
shorter than the gap opening timescale to avoid com¬ 
plicated gap structures (e.g. vortices at the gap edges) 
and other longterm effects (e.g. radial buoyancy waves, 
Richert et al. 2015). To verify that the revealed wave 
mechanics still hold in long terms we have run one sim¬ 
ulation for 120 orbits (SM2IS01ong in Table 1), which 
will be discussed in §6.1. A constant a viscosity with 
a = 10“^ has been applied in our main sets of simula¬ 
tions. 

At inner and outer boundaries, all quantities are fixed 
at the initial states. For a numerical code using Go¬ 
dunov scheme which calculates the flux by decomposing 
wave characteristics at the left and right grid cells, such 
boundary can absorb wave characteristics coming from 
the active zones and limit waves traveling from the ghost 
zones to active zones. Such boundary condition shows 
little wave reflection and is similar to the non-reflecting 
boundary condition (Godon 1996) used in the FARGO 
code (Masset 2000). The detailed code comparison is 
given in Appendix G of Zhu et al. 2014. 

3.2. 2-D simulations 

Gompared with the 3-D simulations in the next sub¬ 
section, 2-D simulations allow us to study density wakes 
in a bigger domain using a higher numerical resolution. 
The initial radial profile of the disk is 

I]o(i?,</>) = I]o(i?o)(T^ (10) 

/ D \ 

To(i?,0)=ro(i?o)f;^] . (11) 

We choose Rq = 1, Ilo(^o) = 1, and Tq{Rq) = 0.01 to 
make {H/R)r=r^ = 0.1. 

Gylindrical coordinates have been adopted. To make 
every grid cell have equal length in the radial and az¬ 
imuthal direction throughout the whole domain, the 
grids are uniformly spaced in log(R) from R =0.2 to 
10, and uniformly spaced from 0 to 27r in the <p direc¬ 
tion. Our standard resolution is 1280 in the R direction 
and 2048 in the (j) direction, which is equivalent to 32 
grids per at R = 1 in both directions. In Table 1, 
2-D runs are denoted with a “G” in front of the model 
names, while 3-D runs are denoted with a “S” in front of 
the model names. 


3.3. 3-D simulations 

To study the 3-D structure of density wakes/shocks, 
we have run 3-D hydro dynamical simulations in spher¬ 
ical polar coordinates except for one case in cylindrical 
coordinates. The initial density profile of the disk at the 
disk midplane is 

po{R,z = Q)= po{Ro,z = Q)(^^^ , (12) 

and the temperature is constant on cylinders 

To(i?,^)=To(i?o)(Ty . (13) 


We want to emphasize that R should not be confused 
with r. In this paper, we use (R, 0, z) to represent po¬ 
sitions in cylindrical coordinates while using (r, 6>, (j)) for 
spherical polar coordinates. <p represents the azimuthal 
direction (the direction of disk rotation) in both coordi¬ 
nate systems. Gonsidering the disk structure is more nat¬ 
ural to be described in cylindrical coordinates, we have 
transformed 3-D simulation results from spherical polar 
coordinates to cylindrical coordinates. Most results pre¬ 
sented below are plotted in cylindrical coordinates with 
R representing the distance to the axis of the disk, even 
though most simulations are carried out in spherical po¬ 
lar coordinates. 

Hydrostatic equilibrium in the r — 0 plane requires that 
(e.g. Nelson et al. 2013) 


po{R,z) = po{Ro,z = 0) 



and 


exp 


GM 




- 

(14) 


1 

R 


^4(R, z) = Dr 


(p + 9) ( § 


+ (1 + 9)- 


qR 




1/2 


(15) 


where = ^JpJp is the isothermal sound speed at R, 

Dr = ^JGM^/R^ ^ and H = Cq/Dr as defined before. 

We choose p = —2.25 and q = —1/2 in our main sets 
of simulations so that Sq oc similar to 2-D sim¬ 

ulations above. po{Rq,z = 0) is 1 and H/R is 0.1 at 
R = Rq. The grids are uniformly spaced in log(r), 6>, <p 
with 256x128x688 grid cells in the domain of [log(0.3), 
log(3)] X [7r/2-0.6, 7r/2+0.6 ]x[0, 27r] for the main sets of 
simulations. In runs with Tcooi = 1 and 100, the cool¬ 
ing time decreases exponentially beyond z = 3i^ with 
Tcooii^) = ^coozexp(—— 3^)) to mimic fast cool¬ 
ing at the disk surface (D’Alessio et al. 1998). Numer¬ 
ically, this treatment also maintains better hydrostatic 
equilibrium at the disk surface. 

The boundary condition in the 0 direction is chosen 
that Vr = Vo = 0 in the ghost zones. We set and T in 
the ghost zones having the same values as the last active 
zones. Density in the ghost zones is set to be 


p{Sg) = p{^a) 


sin(6>g) 

sin(6»a) 


vl/T 


(16) 


to maintain hydrostatic equilibrium in the 0 direction, 
where Og and 9a are the 0 coordinates of the ghost and 
last active zones. We have also tried the boundary con¬ 
dition which sets the quantities in the ghost zones as the 
initial values, and found that the results are not affected 
by the choice of boundary conditions. 

To further study the numerical effect from the bound¬ 
ary, we have applied a wave damping zone (de Val-Borro 
et al. 2006) operating at both r and 0 boundaries in 
the run SM2IS01ong. The wave damping zone in the r 
direction is from Rin to 1.2bRin and from O.SARout to 
Rout‘ The damping zone in the 9 direction is from 9in 
to 9in ^0.1 and from 9out - 0.1 to And Rin, Rout, 
9in, and 9out are the boundary of the simulation domain. 
In these damping zones, the physical quantities are re¬ 
laxed to the initial states on a timescale varying from 
infinity at the damping zone edge 1.2bRin, t).MRout, 
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Oin + 0 . 1 , and Oout — 0.1 to a timescale of 0.1 orbit at 
the boundary of the simulation domain. The damping 
zone gradually damps waves traveling to the boundary 
of the simulation domain. Besides the wave damping 
zone and the long timescale, the run SM2IS01ong is dif¬ 
ferent from our main set of simulations in other ways. 
We ramp the planet mass linearly for 10 orbits in run 
SM2IS01ong to test how our results will be affected if we 
insert the planet gradually in the disk. We also choose 
a = 0 in run SM2IS01ong to confirm that the small vis¬ 
cosity {a = 10 “"^) in the main set simulations will not 
affect the results. 

To verify our results, especially for the 3-D structure 
of spiral shocks, we have used Athena with cylindrical 
coordinates (run SM 2 ISOcyl) to carry out a 3-D simula¬ 
tion having the same disk set-up as in SM2ISO. This 
Athena simulation (SM 2 ISOcyl) is different from the 
Athena++ simulation (SM 2 ISO) in several ways. First, 
SM 2 ISOcyl uses the Corner Transport Upwind Integra¬ 
tor which is different from the Van-Leer Integrator used 
in Athena++. Second, cylindrical coordinates with uni¬ 
form radial grids have been used in SM 2 ISOcyl while 
spherical-polar coordinates with logarithmic radial grids 
have been used in SM2ISO. The detailed disk set-up 
and the boundary conditions can be found in Zhu et 
al. (2014). The results will be presented in §6.1. Over¬ 
all, SM 2 ISOcyl confirms the results in SM 2 ISO, which 
greatly limits the chance that our results are due to nu¬ 
merical artifacts. 

4. THE SHAPE OF SPIRAL WAKES 

Volume rendering of simulation SMIISO is 

shown in Figure Jp is the density difference between 
10 orbits and the initial condition, and po is the initial 
density at that position. Thus, 5p /highlights the den¬ 
sity perturbation (e.g. spiral shocks) in the disk. In this 
paper, we use “spiral shocks” to refer to peaks of the 
density wakes and are associated with spiral arms seen 
in observations. It is apparent in Figurethat the spiral 
shocks are not perpendicular to the disk midplane and 
they have complicated 3-D structure. In the figure, both 
the inner arms inside the planet and the outer arms out¬ 
side the planet curl towards the central star at higher 
altitudes. This curled 3-D shock structure will be stud¬ 
ied in more detail in §5, while in this section we focus on 
the shape of the spiral wakes in the horizontal plane. 

Figure shows the shape of the spiral wakes in both 
2-D and 3-D simulations. The x-axis is plotted in log i?, 
so that the pitch angle of the spiral wake can be easily 
estimated by using its slope in the figure (d log R/d (j) 
=tan /3). 

When a very low mass planet is present in a disk, it 
excites density waves that are in the linear regime. The 
linear theo^ for density waves in the 2-D R — cj) plane 
(Equation can accurately describe the shape of the 
excited spiral wakes in 2-D simulations. This is demon¬ 
strated in the ujmer left panel (CMIISO) of Figure 
where Equation ® with a = 3/2 and 77 = 1/4 fits the 
peak of the density wakes very welQ Even at the mid- 

® Strictly speaking, even with Mp = 0.01 Mj the excited density 
wakes become weak shocks at R = 0.4 and 1.6 according to Equa¬ 
tion But the shocks are very weak and do not move away from 
the trajectory predicted by Equation § significantly. 


plane of 3-D simulations. Equation ([^ still provides a 
good fit to the density wakes (the upper middle panel). 

However, the shape of the spiral arms at the disk sur¬ 
face is affected by the 3-D structure of density wakes. 
At the disk surface in 3-D simulations (even in the linear 
regime, shown in the upper right panel of Eigurej^, both 
inner and outer arms are at smaller R than Equation © 
due to the tilted shock shape in Eigure When these 
shocks are far away from the planet, they are more tilted 
towards the central star at higher altitudes. This leads to 
the inner spiral arms becoming slightly more open (with 
a larger pitch angle) and the outer spiral arms becom¬ 
ing slightly less open (with a smaller pitch angle) than 
Equation Q would predict. 

When the planet has a mass larger than Mth (mid¬ 
dle and bottom panels of Eigure [^, it can launch spiral 
shocks immediately around the planet, and the shape of 
spiral shocks can deviate from the trajectory predicted 
by Equation ([^ significantly. Shocks excited by a more 
massive planet deviate from linear theory more and they 
have larger pitch angles. As shown in Figure [^, spiral 
shocks in the 6 Mj cases (bottom panels) are more open 
and deviate from the prediction of linear theory (dotted 
curves) more than shocks in the 1 Mj cases (middle pan¬ 
els). The deviation from the linear theory has also been 
seen in previous simulations, e.g.. Figure 2 and 10 of de 
Val-Borro et al. (2006), but it has not been explored 
and the physical reason for the deviation is left to be 
unexplained. 

This deviation from linear theory shown in Figure 
is consistent with the predictions from the weakly non¬ 
linear density wave theory by Goodman & Rafikov (2001) 
and Rafikov ( 2002 ). In weakly non-linear theory, the spi¬ 
ral shock can emand in both azimuthal directions away 
from Equation Q, and, at each radius, the shock density 
profile along the azimuthal direction is N-shaped (Figure 
2 in Goodman & Rafikov 2001). The N-shaped shock 
profile expands in the azimuthal direction at a speed 
which is proportional to the normalized amplitude of the 
shock ((Es/ioc/c — So)/^o)- Thus, the higher is the planet 
mass, the stronger are the shocks and these shocks ex¬ 
pand faster away from Equation §. Then, the spiral 
shock has a larger pitch angle as a result. 

Similar to the 0.01 Mj case, the inner spiral shocks 
in 1 and 6 Mj cases are even more open at z = 3H 
than at the midplane, while the outer arms become less 
open at the disk surface. For outer spiral arms, this 3-D 
effect compensates the increased pitch angel due to the 
shock expansion, and coincidently the outer arms almost 
overlap with the prediction from linear theory. 

Another important feature shown in Figure is that, 
besides the primary inner arm which originates from the 
planet, a secondary inner spiral arm appears with some 
azimuthal shift from the primary arm. For some cases 
(e.g. CM 2 ISO), we even see a tertiary arm at the very 
inner disk. The secondary spiral arm has also been seen 
in previous simulations having massive planets, e.g. Fig¬ 
ure 2 in Kley (1999) and Figure 10 of de Val-Borro et al. 
(2006). However, it has hardly been explored in earlier 
simulations. Figure shows that, even with a very low 
mass planet (0.01 Mj, the upper panels of Figure!^, an¬ 
other density peak (the secondary arm) emerges close to 
the primary inner arm with low density region (the rar¬ 
efaction wave of the primary arm) in between. After the 


Q,1 0.2 




Fig. 1. — Volume rendering of dpjpQ for SMIISO. The disk has been sliced through the midplane and meridian plane to show the 3-D 
shock structure. Spiral shocks have been excited by the planet, and the shocks curl towards the central star at the disk surface. 
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Fig. 2.— SpjpQ for CMIISO, CM2ISO, CM3ISO (left panels), and SMIISO, SM2ISO, SM3ISO at the disk midplane (middle panels) and 
2 ; = (right panels). To make the spiral shocks stand out in each panel, we have scaled Spj pQ to a reference value (((5p/po)r ) in each 
panel. ((5p/po)r are 0.01, 0.01, and 0.02 from left to the right panel in the first row, 0.3, 0^ and 0.6 in the second row, and 1, 1, and 2 in 
the third row. The black dotted curves are the spiral wakes from linear theory (Equation]^ with o; = 3/2 and p = 1/4). When the planet 
is more massive, the spiral shocks have larger deviations from the prediction of linear theory. Due to the 3-D structure of the shocks, the 
inner spiral shocks become more open and the outer shocks become less open at 2 ; = 3H compared with the shocks at the midplane. The 
color bar is uniform but it has different scale in each plot to highlight the shock structure. 

secondary inner arm appears close to the primary arm, it can become shock during the propagation and later 
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it will become N-shaped which is similar to the primary 
arm. Then in some cases, a tertiary inner arm appears 
at the rarefaction wave part of the secondary arm. For 
0.01 Mj cases, the primary and secondary inner arms are 
separated hy Sc/) ^ 1 at R = 0.3. When the planet gets 
more massive, the secondary inner arm is excited earlier 
and the separation between the primary and secondary 
arm increases. In 1 Mj cases, the two inner arms are 
roughly separated by (50 ^ 2, and in 6 Mj cases, the 
two arms are roughly separated by (50 ^ 3. This has 
an important application that we can use the separation 
between two arms to estimate the mass of the embedded 
planet. Similar to the primary inner arm, the secondary 
inner arm is also stronger at the disk surface than at the 
disk midplane. For outer arms, the secondary arm ap¬ 
pears in disks that have massive planets embedded, but 
the secondary outer arm is less apparent than the pri¬ 
mary outer arm. More discussions on the 3-D structure 
of secondary arms will be presented in §5. 

Spiral wakes/shocks are slightly more open in a disk 
whose EoS is not isothermal (Figure]^. This is because 
density waves propagate slightly faster in a fluid with a 
non-isothermal EoS than in a fluid at the same tempera¬ 
ture with the isothermal EoS (e.g. Goodman & Raflkov 
2001). In Eigure[^ even with a moderate cooling rate 
{Tcooi = 1): Ihe spiral wakes excited by a low mass planet 
(0.01 Mj) can only be fitted by Equation Q using a 
larger disk scale height that is calculated witn the adia¬ 
batic sound speed instead of the isothermal sound speed 
{hp in Equation]^ is thus Cg^adi/R^ = \/lCs^iso/R^)' H 
is a little bit surprising that the spiral shape in disks with 
Tcooi = 1 follows the spiral shape in adiabatic disks. We 
think this is due to the short timescale for the flow to 
move across the spiral wake. The spiral wake has a typ¬ 
ical width smaller than the disk scale height, while the 
background flow moves across the wake at nearly Kep- 
lerian speed. Thus, when the fluid travels in and out of 
the wave/shock, its response time is much smaller than 
the orbital time and thus behave adiabatically. 

Other aspects of the spiral shocks in non-isothermal 
cases are similar to the isothermal cases, e.g., a higher 
mass planet excites a more open spiral shock, and the 
inner spiral shocks become slightly more open at higher 
altitudes. 

To illustrate the shape of the spiral shocks in the phys¬ 
ical space, we plot the relative density perturbation in 
Cartesian coordinates in Eigure Clearly, the more 
massive the planet is, the more the spiral shocks devi¬ 
ate from linear theory. Two well separated inner arms 
are also apparent when the planet mass is large, and the 
separation between these two arms is larger when the 
planet is more massive (comparing the middle and bot¬ 
tom panels in Eigure®. 

To see how successml the weakly non-linear density 
wave theory of Goodman & Raflkov (2001) and Raflkov 
(2002) fits the shape of the shocks in simulations, we plot 
in Eigure the density contour and density profiles along 
the azimuthal direction atR = 0.3, 0.5, 1.5 and 2 for run 
CM2ISO. In the density contour panel, we can clearly see 
that the secondary arm appears at the edge of the low 
density rarefaction wave region of the primary arm. Af¬ 
ter the secondary arm propagates for a short distance, a 
tertiary arm appears close to the low density rarefaction 
wave region of the secondary arm. In the panels show¬ 


ing the density profiles, we follow Goodman & Raflkov 
(2001) and Raflkov (2002) to shift the density profiles so 
that (/)' = 0 corresponds to the wake position from linear 
theory (black curve in Eigure]^ or Equation [^. The den¬ 
sity profiles clearly show that the shocks are N-shaped. 
A rarefaction wave follows the shock front, and Sp at the 
rarefaction wave region can be negative before it merges 
to the background flow. The shock fronts deviate from 
0' = 0, and due to the shock expansion, the deviation is 
larger when the shock is further away from the planet. 
Under the shearing-sheet approximation, the amplitude 
and width of the N-shaped shock scale as |R — Rp\~^^^ 
and \R — Rp\^^^ at \R — Rp\ ^0 based on the weakly 
non-linear density wave theory of Goodman & Raflkov 
(2001). In a global disk spanning a large range of radii, 
the amplitude and width of the N-shaped shock scale as 
and where t is given in Equation 43 of Raflkov 
(2002). Eor our disk parameters, we have 


t (X 



^3/2_i|3/2^-13/8^^ 


(17) 


Since the middle point of the N-shaped shock is at 0' ^0 
(Raflkov 2002), we expect that in our global simulations 
the azimuthal distance between the shock front and the 
path predicted by Equation ® (0'=O) should also scales 
as 

To test this prediction, we have measured the shock 
positions at R = 0.5 and 1.5 in Eigure which are 
0' = —0.9 and 0.43 respectively. These two positions 
are labeled as the dashed lines in R = 0.5 and 1.5 pan¬ 
els. Then we calculate the shock positions at R = 0.3 
and 2 to be -1.7 and 0.9, using their positions at R = 0.5 
and 1.5 together with the scaling relationship where 
t at R=0.3, 0.5, 1.5, and 2 are calculated from Equa¬ 
tion |T^ These predicted shock positions are labeled as 
the dashed lines in R = 0.3 and 2 panels. We can see 
that they agree with the actual shock positions in the 
simulation very well. This confirms that the shock posi¬ 
tions are determined by the non-linear expansion of spiral 
shocks. Using the same approach, we have calculated the 
non-linear shock position at every radius for R < 1 and 
R > 1 with the normalization based on shock positions 
at R = 0.5 and 1.5. This new predicted shock shape is 
plotted as the dotted curve in the left panel of Eigure 
Despite some offset at R close to the planet, which is 
expected since the simple relationship holds only 

when \R — Rp\ ^0, a good agreement has been achieved 
between the non-linear density wave theory and the sim¬ 
ulations. Note that in this comparison, we did not cal¬ 
culate the shock strength directly from non-linear theory 
and compare its amplitude with simulations, instead we 
use the scaling relationship to verify the propagation of 
the shock. In future, direct comparison is desired when 
we have a more complete non-linear theory which can 
calculate the shock excitation directly. 

Although the primary arm can be fitted by the weakly 
non-linear density wave theory, the excitation of the sec¬ 
ondary (or even tertiary) arm still lacks a good theoret¬ 
ical explanation. It may be related to the low m mode 
(e.g. m = 2, 3, similar to disks in binary systems) or 
some non-linear wave coupling. Eigure [^suggests that a 
secondary spiral arm is excited at the other end of the 
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Fig. 3.— The same as Figure|^but for CMITI, CM2T1, CM3T1 (left panels), and SMITI, SM2T1, SM3T1. {5p/po)r are 0.005, 0.005, 
and 0.01 from left to the right panel in the first row, 0.3, 0.3, and 0.6 in the second row, and 0.8, 0.8, and 2.4 in the third row. The black 
squared dots represent the linear theory using isothermal sound speed while the black plus sign dots use the adiabatic sound speed. 




Fig. 4.— The same as Figurej^but in Cartesian coordinates. Please note that ((5p/po)r is different in each panel to make the spiral arms 
stand out. 

N-shaped primary shock (the R = 0.5 panel). After it 


is excited, it steepens to shocks and becomes another N- 
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Fig. 5.— ^p!Pq for CM2ISO (left color panel) and density cuts across the 0 direction at R =0.5, 0.3, 1.5, and 2.0. The solid black curve 
in the left color panel labels the shock position from linear theory (Equation!^. In the density cut plots, the density profile has been 
shifted so that the shock position from linear theory is at (j)' = 0. The dashed nnes in i?=0.5 and 1.5 plots label the shock fronts, while 
the dashed lines in i^=0.3 and 2 plots label the predicted shock fronts from weakly nonlinear theory. In the left color contour panel, the 
predicted shock front from nonlinear theory is labeled as the dotted curve. 


shaped shock later (the R = 0.3 panel). Its shock front 
can even travel into the rarefaction wave of the primary 
arm (e.g. in the R = 0.3 panel, the secondary shock 
is almost at 0'=O where the rarefaction wave of the pri¬ 
mary arm should reside.). Unlike the primary arm which 
already dissipates when it travels inward from = 0.5 
to 0.3, this secondary arm is excited later and becomes 
stronger from i? = 0.5to0.3. Ati? = 0.3, the secondary 
arm is even stronger than the primary arm. By compar¬ 
ing = 0.5 and = 0.3 panels, we also notice that the 
secondary arm almost keeps the same azimuthal separa¬ 
tion with the primary arm (A0' ^1.7) during its prop¬ 
agation. Furthermore, at i? = 0.3, a tertiary arm starts 
to appear at the other end of the N-shaped secondary 
arm. 

To summarize the results in this section. Figure 
shows the pitch angle from the linear theory and those 
measured in numerical simulations. When the planet 
mass is low (e.g. CMIISO, the green dots), the mea¬ 
sured pitch angle agrees with the linear theory. When 
the planet mass increases, the pitch angle also increases. 
For the 6 Mj case, the measured pitch angle of the spi¬ 
ral wake in the hp = 0.1 disk is close to the pitch angle 
predicted in a much thicker disk {hp = 0.2) using the 
linear theory. The secondary and even the tertiary arms 


have similar pitch angles as the primary arms. Thus, if 
we know the disk thermal structure very well, we can use 
the deviation of the measured pitch angle from the linear 
theory to estimate the embedded planet mass. 

5. THE 3-D STRUCTURE OF SPIRAL WAKES 

Since near-IR scattered light observations only probe 
the shape of the disk surface, 3-D structure of spiral 
shocks can affect the observational signatures of these 
spiral shocks. Intuitively, we would expect that the spiral 
shocks have complicated 3-D structure. First, the wave 
excitation must have 3-D structure since, at the same 
R in the disk, the distance between the planet and the 
disk surface is larger than the distance at the midplane, 
and the force is thus weaker at the disk surface. Second, 
the wave propagation may have 3-D structure consider¬ 
ing the disk becomes thinner at smaller R. Waves/shocks 
are more converged when they propagate inwards. They 
can also channel to the disk surface (Lubow & Ogilvie 
1998, Bate et al. 2002, Lee & Gu 2015), and, during 
their propagation from the high density region (e.g., the 
midplane) to the low density region (e.g., the disk sur¬ 
face), the amplitudes of perturbations have to increase 
to conserve the wave action. Since the amplitudes of 
perturbations determine when the waves will break into 
shocks, the dissipation can also be quite different between 
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Fig. 6 .— The pitch angle based on the linear theory (Equation 
with OL = 3/2 and 77 = 1/4) assuming hp=0.1 (solid curve), 
U.2 (dotted curve), and 0.05 (dashed curve) compared with those 
measured in numerical simulations. The red, blue, and green dots 
are measured from CM3ISO, CM2ISO, and CMIISO respectively. 
The circles, triangles, and crosses are the pitch angle of the primary, 
secondary and tertiary arms. 
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Fig. 7. — Waves excited on a 3-D disk by a low mass planet 
in run SMIISO (left panels) and STHIN (right panels). Density 
perturbations with n=0,2, and 4 Hermite components for the m=10 
Fourier mode are displayed. The real (imaginary) part of 77 is shown 
with a dotted (solid) curve. Y-scales are different for different n 
modes. 

the surface and the midplane. All these effects can con¬ 
tribute to the 3-D structure of spiral waves/shocks. 

Due to these complicated effects, it is difficult to de¬ 
velop an analytic theory to study the planet-induced 3-D 
shock structure, and we rely on numerical simulations to 
study such structure. However, before delving into the 
highly nonlinear shock regime, we can use the linear the¬ 
ory developed in Tanaka, Takeuchi & Ward (2002) to 
estimate the 3-D effect of the density waves. Following 
their theory for locally isothermal disks, the structure of 


the waves in the 2 ; direction can be studied with Her¬ 
mite polynomials (i^^(Z)). We first expand perturbed 
quantities {r]) from our simulations into Fourier series 


T] = y^Re 


rime 




(18) 


where the Fourier components 77 ^ are complex functions 
of R and 2 ). Then, 77 ^ can be further expanded with 
Hermite polynomials in the 2 : direction. 


— ^ ^ ^ ^ Re Tjrfji nHYii^Z^e 


m=0n=0 




(19) 


where Z is the normalized height as Z = zlH{R)^ and 
the first three Hermite polynomials are 

Ho{Z) = l, H^{Z) = Z, H2{Z) = Z^-1. (20) 

By using the normal orthogonal relation between we 
have 


Vm,n 



^"/^HniZ)7l^dZ. 


( 21 ) 


We can use r]rn,n at different n to estimate the relative 
importance of different Hermite components. To com¬ 
pare with Figure 2 in Tanaka et al. (2002), we show 
m = 10, n = 0,2,4 Fourier-Hermite components for the 
perturbed density {6p/po^m FigurelH With the same 
parameters, the right panS^of Figure [^is very similar to 
Figure 2 in Tanaka et al. (2012). By comparing the left 
and right panels in Figure we can see that although 
higher-order Hermite components have similar strength 
between thin and thick disks, n=0 Hermite component is 
much weaker in a thick disk, suggesting that the wakes 
in a thick disk have more significant 3-D structures than 
wakes in a thin disk. 

Figure suggests that higher-order vertical compo¬ 
nents can dominate the disk structure at the atmosphere. 
Although it shows that the n = 4 component ( 7710 , 4 ) is 10 
times weaker than the n = 2 component ( 7710 , 2 ), and the 
n = 2 component is 10 times weaker than the 77 . = 0 com¬ 
ponent (which led Tanaka et al. 2002 to conclude that 
most of the angular momentum excited by the planet will 
be carried by two dimensional free waves), the base func¬ 
tion (Hermite polynomials) at ^ = 3i^ gets ^ 10 times 
larger sequentially from Hq to H 2 and to Thus, 

^ 10 , 2^2 and 7710 , 4 i ^4 are still comparable with 7710 , 0 ^ 0 - 
The density structure at the disk atmosphere can be sig¬ 
nificantly affected by higher-order vertical modes. 

Although the modal analysis is useful to verify numeri¬ 
cal simulations and can be suggestive on the relative am¬ 
plitudes of various modes, it is the 3-D structure in the 
real space that determines the observational signatures 
of waves/shocks. 

By studying the shock structure in real space, we first 
find that the 3-D shock structure is dramatically differ¬ 
ent between inner and outer arms. For the inner arms, 
the density perturbation is much larger at the disk sur¬ 
face than at the disk midplane. The 3-D structure of 
the inner spiral arms at R = 0.5 is shown in Figure 


^ Since the disk is isothermal locally, the density perturbation 
is quite similar to the enthalpy perturbation, and can be compared 
with Figure 2 in Tanaka et al. (2012). 
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Fig. 8 .— At the disk radius of R = 0.5, density profiles along 
the azimuthal direction at the disk midplane (dotted curves) and 
2 ; = 3H (solid curves). {p)cf) is the density averaged over the az¬ 
imuthal direction. Simulations with different planet masses (dif¬ 
ferent columns) and equations of state (upper panels: isothermal, 
bottom panels: adiabatic with Tcool = 1) have been shown. 


At R = 0.5 and 2 : = 3H (solid curves), the differences 
between the maximum and minimum density in the loga¬ 
rithmic scale are 0.015, 0.4, and 1.3 for Ml, M2, and M3 
cases respectively, in comparison with 0.004, 0.2, and 0.4 
at the disk midplane (dotted curves). At the same radius 
{R = 0.5), the position of the wakes in non-isothermal 
disks (lower panels) are at smaller (p — (pp compared with 
those in isothermal disks (upper panels). This is because 
the wakes are more open in non-isothermal disks, as dis¬ 
cussed in §4, . 
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Fig. 9.— Vr (left panels) and vq (right panels) at 2 : = IH (dot¬ 
ted curves) and 2H (solid curves) for SMIISO (upper panels) and 
SMITI (bottom panels), vq is positive when the motion is towards 
the disk surface. 


The secondary inner spiral arms/shocks are also more 
prominent at the disk surface than at the disk midplane. 
At the disk midplane, the secondary arms have lower am¬ 
plitudes compared with the primary arms (dotted curves 
in Figure]^, while at 2 ; = 3i7, the secondary arms have 
almost the same amplitudes as the primary arms (solid 
curves). The large amplitude of primary and secondary 
inner arms at the disk surface is due to the corrugated 


motion in the vq direction. In Figure which is in the 
corotating frame with the planet, the disk material ffows 
in the direction from the left side to the right side of 
the figure. Before meeting with the shock, the disk is 
in vertical hydrostatic equilibrium with the background 
density and vq = t) (Figure [^. After the shock, the 
disk material loses angular momentum and moves in¬ 
wards with Vr < ^ (Figure [^. At the same time, vq 
also becomes negative, compressing the disk material at 
the midplane. This downward motion decreases the den¬ 
sity of the rarefaction wave at ^ = 3H. Before meeting 
the secondary shock, vq starts to increase and becomes 
positive, leading to a higher density at the disk surface. 
At the secondary shock, vq reaches the maximum posi¬ 
tive velocity and leads to the highest density at the disk 
surface for the secondary shock in Figure This corru¬ 
gated motion, first negative and then positive vq^ leads 
to an enhanced contrast between the spiral shock and 
the rarefaction wave after the shock. 




Fig. 10. — Similar to Figurebut at = 2. 


On the other hand, the density perturbation of outer 
spiral arms is similar between the disk surface and the 
disk midplane, especially for isothermal disks, as shown 
in Figure At = 2, regardless of height, the dif¬ 
ferences between the maximum and minimum density 
in the logarithmic scale are both 0.002 for SMIISO, 0.1 
for SM2ISO , and 0.3 for SM3ISO. This lack of vertical 
variation is also reflected in Figure [TT] where vq is very 
small compared with Vr {vq is almost two orders of mag¬ 
nitude smaller than Vr). Thus, the density structure of 
the outer spiral is mainly determined by Vr and in the 
horizontal plane in stead of the corrugated motion in the 
0 direction. 


For non-isothermal runs (bottom panels in Figure 10), 
the density perturbation of the outer spiral arms at the 
disk surface is slightly higher than the perturbation at 
the disk midplane. Disk material flows from the right 
hand side of the figure to the left hand side in Figure pT| 
and When it meets the shock, it develops a vq to- 
wardsthe disk surface, which enhances the density at the 
disk surface. Although it is tempting to contribute such 
difference between isothermal and non-isothermal runs 
to the nonlinear hydraulic jumps (shock bores) (Boley & 
Durisen 2006), such disk structure also appears even in 
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Fig. 11.— Similar to Figure but at R = 2. 

the linear regime for the 0.01 Mj case, implying that it 
may be a linear effect and related to the eigenfunctions 
of the 3-D waves excited by the planet. 

Finally, to illustrate the increase of the density per¬ 
turbation with height in disks and the qualitative dif¬ 
ference between inner and outer arms, we plot the 
relative density perturbation along the radius at dif¬ 


ferent heights {z =0, 1, 2, 3, 4 H) in Figures 12 
and The relative density perturbation is delink 

as pmax{R,z)lpmin{R,z) - 1, where pmax{R,z) and 
Pmin{R^ are the maximum and minimum density along 
the azimuthal direction (0 = [0,27r]) at the fixed R and 
2 ;. We can see that the inner arms and outer arms are 
qualitatively different. For inner arms, the relative den¬ 
sity perturbation is getting larger at higher altitudes, 
and the perturbation can increase by more than a fac¬ 
tor of 10 from the midplane to the disk surface. For the 
outer arms, the relative density perturbation is almost 
unchanged between the midplane and the disk surface in 
isothermal disks (Figure and only increases slightly 
from the mid plan e to the oisk surface in non-isothermal 
disks (Figure]!^. Overall, for the inner arms, the large 
density pertuilDation at the disk surface has a significant 
effect on the near-IR observations as shown below. 


6. DISCUSSION 

6.1. Numerics, Longterm Evolution, and Different Disk 

Structures 

As shown in Figure[^ the 3-D spiral shocks are not per¬ 
pendicular to the disk midplane. They curl in a way that 
they are almost along the 0 direction in the spherical- 
polar grid. To quantify the curl, the upper panel of Fig¬ 
ure shows the relative density perturbation at the 
disk midplane (black curves) and at 0 = 7r/2 + 0.35 
(cyan curves) for run SM2ISO (solid curves). Please note 


that the x-axis in Figure 14 is r (the radial position in 


spherical-polar coordinates) instead of R (the radial po 
sition in cylindrical coordinates). The fact that the shock 
density peaks at almost the same r at either the disk mid¬ 
plane OT 0 = 7r/2 + 0.35 demonstrates that the shocks al¬ 
most curl along the 0 direction. Such coincidence makes 
us suspect that they could be numerical artifacts due to 
the adopted spherical-polar grid structure. Thus, we use 
Athena to carry out a similar simulation but using cylin¬ 
drical coordinates (run SM2ISOcyl introduced in §3.3). 


After transforming simulation outputs from cylindrical 
coordinates to spherical-polar coordinates, the relative 
density perturbation for SM2ISOcyl is shown in the up¬ 
per panel of Figure as the dotted curves. The shocks 
are at similar positions as those in SM2ISO, which con¬ 
firms that the curled shock structure is real instead of 
numerical artifacts. 

The physical reason for the curled shock structure is 
unclear. One would suspect that it may be related to 
the radial temperature gradient in the local isothermal 
disk. In such disks, the disk rotates at different angular 
velocities at different disk heights. Such vertical shear 
can change disk dynamics, such as leading to the vertical 
shear instability (Nelson et al. 2013) and it may also lead 
to curled shock structures. To test this idea, we have run 
one simulation with a constant temperature in the whole 
disk (SM2ISOconstT). In such a disk, the disk rotates 
at a constant angular velocity at a given R independent 
on the disk height. The relative density perturbation is 
shown in the bottom panel of Figure [14| However, the 
shock density still peaks at almost the same r at different 
0, implying that the curl of the shock is not caused by 
the vertical shear in disks. Thus, the physical mechanism 
for the curled shock remains unclear and still deserves 
further explore. 

To demonstrate that our derived shock structure is in¬ 
dependent on specific numerical choices in our model, 
we show the density perturbation for run SM2IS01ong in 
Figure p+1 Run SM2IS01ong is different from SM2ISO in 
several aspects (as introduced in §3.3) : 1) the simulation 
runs for 120 orbits instead of 10 orbits; 2) it has a wave 
damping region, 3) the planet mass is ramped up slowly 
over 10 orbits, and 4) it uses zero viscosity. Despite these 
differences, its density perturbation at 15 orbits (bottom 
panel of Figure 15) is quite similar to the density pertur¬ 
bation in SM2ISD (middle panel of Figure [T^. When 
a gap is induced at 120 orbits, the spiral shocks become 
weaker at the disk region close to the gap, while it still 
maintains the same strength at the region far away from 
the gap. We notice that vortices start to develop at the 
outer gap edge {R ^ 1.4) at 120 orbits, which slightly 
enhances the density perturbation at the outer gap edge. 

We have also studied how the shock structure can be 
affected by different disk structures. Figureshows the 
relative density perturbation of the wakes in a thinner 
disk (STHIN). Compared with Figurewe can see that 
the density perturbation increases by a factor of 2 from 
the midplane to 3 H, and a factor of 4 from the midplane 
to 4 in the thin disk, compared with a factor of 4 
and 10 respectively in the thick disk. This suggests that 
the wakes/shocks have more significant 3-D structure in 
a thicker disk. Although this does not favor detecting 
spiral arms in thinner disks in future, we need to keep in 
mind that the thermal mass (Equation is smaller in 
a thinner disk so that the same mass planet corresponds 
to a more massive planet in the scale of thermal mass 
and it will excite stronger density waves. In this case the 
inner arms may still be observable in a thin disk. 


6.2. Near-IR Images 

To understand how the 3-D structure of density 
waves/shocks affects observations, we post-process our 
hydro dynamical simulations with Monte-Carlo radiative 
transfer calculations to generate near-IR scattered light 
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Fig. 12.— Relative density perturbations at different heights for SMIISO, SM2ISO, SM3ISO. pmax and pmin are the maximum and 
minimum density along the circle in the azimuthal direction at given R and 2 ;. The black (blue, cyan, orange, green) curve is calculated at 
the disk midplane (lif, 2if, 3if, AH). The shaded region represents the region where the density perturbation is determined by the planet 
and the circumplanetary region instead of the spiral shock. 
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Fig. 13.— Similar to Figure [T^ but for SMITI, SM2T1, and SM3T1. 


images (Figures 17, 18 and 19). The details on the 


Monte-Carlo radiative transfer calculations are presented 


in Dong et al. (2014, 2015). To assign physical scales to 
our simulations, we assume that the planet is at 50 AU 













































































14 


Zhu et al. 



0.5 1.0 1.5 2.0 2.5 


r 

Fig. 14. — The relative density perturbation at the disk midplane 
(black curves) and at 6 = 7r/2 + 0.35 (cyan curves) for SM2ISO 
(solid curves in the upper panel), SM2ISOcyl (dotted curves in the 
upper panel), and SM2ISOconstT (solid curves in the lower panel). 



R 


Fig. 16.— Density perturbation for run STHIN. Different colors 
represent density perturbation at different disk heights similar to 
Figure!^ 



R 

Fig. 15.— The disk surface density (upper panels) and den¬ 
sity perturbation (lower panel) at 15 (dotted curves) and 120 
(solid curves) orbits for run SM2IS01ong. Different colors in 
the lower panel represent density perturbation at different disk 
heights similar to Figure The disk surface density perturba¬ 
tion (Emacc/Smin-1) IS very similar to the density perturbation 
{ pmax / Pmin -^) at the disk midplane. 

and the central source is a typical Herbig Ae/Be star (2 
Mq) with a temperature of 10^ K and a radius of 2Rq. 
ISM dust grains have been used and their distribution is 
assumed to follow the gas distribution. Following Dong 
et al. (2015), we assume that the gas disk is 0.02 Mq. 
With the gas to dust mass ratio of 100:1, the total dust 
mass is 2x10“^Mq. Further assuming 10% of dust is in 
the form of ISM dust, the total mass of the ISM dust 
is 2 xlO“^M 0 . In this model, Toomre Q parameter at 
50 AU is ^30, so neglecting disk self-gravity is justi¬ 
fied. In MCRT simulations, photons from the central 
star are absorbed/reemitted or scattered by the dust in 
the surrounding disk. The convolved images in Figures 
pT| p! 8 | and are derived by convolving full resolution 
images with a Gaussian point spread function having a 
full width half maximum (FWHM) of 0.06”. This resolu¬ 
tion is comparable with NIR direct imaging observations 


using Subaru, VLT, and Gemini. In the right two panels 
of Figures and we assume that the object is 140 
pc away, while in Figure and the left two panels of 
Figures p! 8 | and \T^ we assume that the distance is 70 pc 
so that tHe inner arms are shown more clearly. 

Both full intensity and polarized intensity images are 
calculated and shown in Figure The spiral arms are 
evident in these images. When the disk is viewed face- 
on, both full intensity and polarized intensity images are 
almost identical except that the full intensity image is 
almost a factor of 2 brighter than the polarized intensity 
image. The polarization fraction, which is defined as the 
ratio between the polarized intensity and full intensity, 
is almost a constant (^0.45). When the disk is viewed at 
some inclination angle, the polarization fraction is not a 
constant due to the dust forward scattering. The depen¬ 
dence of the polarization fraction on disk inclination can 
be useful to determine the disk inclination angle. 

To highlight the importance of the 3-D wave struc¬ 
ture, we have also computed models by only using the 
disk midplane density from simulations, which is labeled 
as 2D-^3D in Figures andIn these models, we as¬ 
sume that the disk is in verti^ hydrostatic equilibrium 
and puff up the midplane density to higher altitudes as 
P = where h/R = O.IR^-^^ (the same 

as the scale height used in 3-D simulations). 

Figures and show that, in 3D models, the inner 
arms are considerably more prominent than the outer 
arms, and normally a secondary arm can be as bright 
as the primary arm. The shape of the inner arms clearly 
deviate from the prediction of linear theory. As expected 
from the non-linear expansion of spiral shocks, the pitch 
angle of the inner spiral arms in the more massive planet 
case (SM3ISO, Figure [l^ is larger than those in the less 
massive planet case (SM 2 ISO , Figure [T8|. The inner 
arms are also quite sharp, while the outer arms are quite 
broad and sometimes indistinguishable from the back¬ 
ground disk. This difference is partly because the sharp 
shock fronts are facing the star for the inner arms, while 
they are facing away from the star for the outer arms. 
The different geometry at the disk surface can greatly 
affect the intensity of the scattered light images (Takami 
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Fig. 17.— Pull intensity (left panels) and polarized intensity after scaled to (middle panels), and polarization fraction (right panels) 
for SM3ISO (Mp = 6Mj) when the disk is viewed at the 60^ inclination angle (upper panels) and face on (bottom panels). 


et al. 2014). We calculate an approximate scattering 
surface, defined as the disk surface where the column 
density is 0.01 (in code units), for the SM3ISO model at 
(j) — (pp = 78^, shown in Figure]^ Clearly, for the inner 
arms, the shock fronts are facing the star, while, for the 
outer arms, the smooth rarefaction waves are facing the 
star. Since the rarefaction waves change gradually with 
radius, they are illuminated by the star more uniformly 
than the shock. Thus, the outer arms appear quite broad. 
However, when the planet mass is not very high (1 Mj 
case), the width of rarefaction waves in the radial di¬ 
rection can be smaller than the size of the observational 
beam, and we won’t be able to distinguish the inner and 
outer arms based on the sharpness of the arms. 

In 3-D models of Figure and the secondary inner 
arms are as apparent as the primary arms, even though 
the primary arms have higher surface density than sec¬ 
ondary arms. This is due to the corrugated motion dis¬ 
cussed above in Figure and which increases the den¬ 
sity of the secondary arms at the disk surface. The sec¬ 
ondary arm is offset from the primary arm with some az¬ 
imuthal angle, as also shown in the surface density plot 
(Figure 1^. This offset is smaller in SMIISO (Figure]!^ 


than that in SM6ISO (Figure I^. The two spiral arms 


are ^ 100*^ apart in SMIISO (Figure 18), while almost 


180^ apart in the more massive planet case (SM6ISO 
Figure 


19). 


By comparing 3D models with 2D^3D models in Fig¬ 
ure and we find that the inner spiral arms are 
more prominent in 3D models, as expected due to inner 
arms’ higher densit y p erturbation at the disk surface in 
3D models (Figure |T^. Even in convolved images, the 
polarized intensity oimner arms is at least twice stronger 
in 3D models than in 2D—)^3D models. 

The outer shocks are not very apparent, and they are 
similar between 3D and 2D^3D models since the density 
perturbation of outer arms is almost height independent 
(§5). As discussed in §4, the outer arms coincidently 


follow linear theory 


is very large as in Figure 


)ry (Equation!^, 
in Figure fl^ th( 


When the planet mass 
e secondary outer arm 


starts to become visible . 

Another noticeable difference between 3D and 2D^3D 
models is that the planet (or the circumplanetary region) 
is bright in 2D^3D models while it is dim in 3D models. 
This is because, when we puff the disk from 2-D to 3-D, 
we have ignored the planet’s gravity so that the higher 
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Fig. 18.— The near-IR polarized intensity maps for inner (left two panels) and outer (right two panels) arms for SM2ISO (Mp = Mj). 
The planet is assumed at 50 AU. In the left two panels, we assume the system is 70 pc away so that the inner arms are shown clearly, while 
in the right two panels we assume the system is 140 pc away. The full resolution images have been convolved with a 0.06” beam to derive 
the convolved images. The upper panels show the images using disk structure directly from 3-D simulations while the bottom panels use 
disk structure assuming the disk is in vertical hydrostatic equilibrium. The dotted curves are the positions of the spiral wake derived from 
linear theory (Equation]^. 
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Fig. 20. — The scattering surface where the column density is 
0.01 (code units) for SM3ISO. Shock fronts are facing the star for 
the inner spiral shocks while the rarefaction waves are facing the 
star for the outer spiral shocks. 

even casts a shadow to the outer disk in 2D->3D mod¬ 
els. In realistic 3-D models, the gravity of the planet has 
been self-consistently included, which pulls the circum- 
planetary material towards the disk midplane and leads 
to a lower density at the disk surface. Thus, the circum- 
planetary region receives less irradiation by the central 
star, and becomes dark. On the other hand, the outer 
disk beyond the planet is better illuminated and thus 
becomes bright instead of being shadowed. 

However, we need to keep in mind that we have ignored 
the luminosity from the planet and the circumplanetary 
disk in our models. Zhu (2015) point out that accreting 
circumplanetary disks can be very bright (^ 0.001 Lq if 
the circumplanetary disk accepts onto Jupiter at a rate 
^ lO“^M 0 /^r). Such high luminosity may be able to 
illuminate the circumplanetary region significantly. We 
may also be able to directly detect such accreting circum¬ 
planetary disks in direct imaging observations operating 
at mid-IR wavelengths. 

6.3. Observational Implications 

Previous works suggest that spiral patterns from di¬ 
rect imaging observations cannot be explained by planet- 
induced spiral arms, since the pitch angle of the spiral 
arm is larger in observations than that predicted by the 
linear spiral wave theory, and the contrast of the spiral 
arm is higher in observations than suggested by the syn¬ 
thetic observations based on two dimensional planet-disk 
simulations. However, most previous works only focus on 
using outer spiral arms beyond the planet to explain ob¬ 
served spiral patterns. As shown in Figure the inner 
arms generally have larger pitch angle than outer arms 
even under the linear density wave theory. Furthermore, 
when the spiral arms become spiral shocks, the pitch an¬ 
gle increases with stronger shocks. We also found that a 
secondary (or even a tertiary) spiral arm, especially for 
the inner arms, is also excited by a massive planet. The 
more massive is the planet, the larger is the separation 
in the azimuthal direction between the primary and sec¬ 
ondary arms. The inner arms also have significant ver¬ 
tical motion, which boosts the density perturbation at 
the disk surface and the intensity contrast in synthetic 
images. 


Thus, we have three independent ways to estimate the 
planet mass based on the spiral patterns in observations: 
1 ) using the deviation between the measured pitch angle 
and the pitch angle predicted in the linear theory (Figure 
[^; 2 ) the existence of the secondary ^iral arm and the 
separation between two arms (Figure!^; 3) the intensity 
contrast of the spiral arms in observation (Figure p! 8 | fc 


19). The first two methods only use the information on 
the shape of the spiral arms and are less likely to be 
affected by the disk vertical temperature structure. In 
a companion paper (Dong et al. 2015), we have shown 
that all these three methods suggest that SAO 206462 
and MWC 758 may harbor massive planets (several to 
tens of Jupiter masses) outside the detected spiral arms. 

The planet-induced spiral arms can efficiently scatter 
light from the central star and are quite bright within 
gaps (Figure fl^. The spiral arms may be responsible 
for some resolved infrared features within the cavity of 
LkCa 15 (Kraus & Ireland 2012 ). 


7. CONCLUSION 

We have carried out two dimensional ( 2 -D) and three 
dimensional (3-D) hydrodynamical simulations to study 
spiral wakes/shocks excited by young planets. Simula¬ 
tions with different planet masses (0.01, 1, and 6 Mj) 
and different equations of state (isothermal and adia¬ 
batic) have been carried out. 


• We find that the linear density wave theory can 
only explain the shape of the spiral wakes excited 
by a very low mass planet (e.g. 0.01 Mj). Spi¬ 
ral shocks excited by high mass planets clearly de¬ 
viate from the prediction of linear theory. For a 
more massive planet, the deviation is more signifi¬ 
cant and the pitch angle of the spiral arms becomes 
larger. This phenomenon can be nicely explained 
by the wake broadening from the non-linear density 
wave theory (Goodman & Rafikov 2001, Rafikov 
2002 ). A more massive planet excites a stronger 
shock which expands more quickly, leading to a 
larger pitch angle. 

• A secondary (or even tertiary) inner spiral arm is 
also excited by the planet. It seems to be excited 
at the edge of the N-shaped primary arm. The 
more massive is the planet, the larger is the sep¬ 
aration between the primary and secondary arm. 
At the disk surface, the secondary inner arm can 
be as strong as the primary arm. The secondary 
inner arm almost keeps the same azimuthal sepa¬ 
ration with the primary arm at every radius in the 
disk. The excitation mechanism for the secondary 
(tertiary) arm deserves further study. 

• The spiral shocks have significant 3-D structure. 
They are not perpendicular to the disk midplane. 
They are curled towards the star at the disk sur¬ 
face. This further increases the pitch angle of the 
inner arms at the disk surface, but reduces the 
pitch angle of the outer arms at the disk surface. 
For outer arms, this effect compensates the in¬ 
creased pitch angle due to wake broadening. Even¬ 
tually at the disk surface, the shape of outer spiral 
arms still roughly follows the prediction of linear 
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theory, while the inner arms are considerably more 
opened than predicted by the linear theory. 

• The inner spiral shocks also have significant vertical 
motion. The corrugated motion increases the den¬ 
sity perturbation of the inner spiral arms by more 
than a factor of 10 at 2 : ^ 3 — AH compared with 
the perturbation at the disk midplane. This can 
dramatically increase the contrast of the spiral pat¬ 
terns in near-IR scattered light images. The outer 
spiral shocks have little vertical motion in isother¬ 
mal disks. With a non-isothermal EoS, there are 
some vertical motions for the outer arms, which 
can make the outer arms more apparent. 

• We have combined our hydro dynamical simulations 
with Monte-Carlo radiative transfer calculations to 
generate near-IR scattered light images. We find 
that the inner spiral arms are prominent features 
that are observable by current near-IR imaging fa¬ 
cilities. Besides the apparent spiral patterns in 
some transitional disks, planet-induced spiral arms 
may also be responsible to some marginally de¬ 
tected near-IR features within cavity of some tran¬ 
sitional disks. We further demonstrate that inner 
spiral arms in synthetic near-IR images using full 
3-D hydro dynamical models are much more promi¬ 
nent than those based on 2-D models assuming 
hydrostatic equilibrium, consistent with the 3-D 
structure of the inner arms. On the other hand, 
the outer shocks are not very apparent and they 
are similar in synthetic images using 3D and 2D 
models. This indicates the need to model obser¬ 
vations (especially for inner arms) with full 3-D 
hydrodynamics. 

• The different geometry between the inner and outer 
arms also affects their appearance in near-IR im¬ 
ages. The sharp shock fronts of the inner arms 
face the central star directly, producing sharp nar¬ 
row spiral features in observations. On the other 
hand, for the outer arms, the smooth rarefaction 


waves face the central star, producing broad and 
dimmer spiral features. 

• In near-IR images, the circumplanetary region is 
very dim since the planetary gravity reduces the 
density at the disk atmosphere. However, the disk 
region behind the planet can be better illuminated 
and becomes bright. 

In the Appendix, we have shown that buoyancy reso¬ 
nances are confirmed in global adiabatic simulations even 
if the disk has a moderate cooling rate. They can lead to 
sharp density ridges around the planet, which may have 
observational signatures. 

Overall, spiral arms (especially inner arms) excited by 
low mass companions are prominent features in near-IR 
scattered light images. Most importantly, we have three 
independent ways to infer the companion’s mass: I) the 
pitch angle of the spiral patterns, 2 ) the separation be¬ 
tween the primary and secondary arms, and 3) the con¬ 
trast of the spiral patterns. 

In a companion paper (Dong et al. 2015), we have com¬ 
bined MCRT and hydro dynamical simulations from this 
paper, and shown that planet-induced inner arms can 
explain spiral patterns revealed by recent near-IR direct 
imaging observations for SAO 206462 and MWC 758. 
We want to caution that our proposed model can not 
explain the gaps at the inner disks. Other mechanisms 
(e.g. another planet or disk photoevaporation) have to 
be invoked to explain the inner gaps. 

All hydro dynamical simulations are carried out using 
computer supported by the Princeton Institute of Com¬ 
putational Science and Engineering, and the Texas Ad¬ 
vanced Computing Center (TACC) at The University of 
Texas at Austin through XSEDE grant TG- ASTI30002. 
This project is supported by NASA through Hubble Eel- 
lowship grants HST-HE-5I333.0I-A (Z.Z.) and HST-HE- 
51320.OI-A (R.D.) awarded by the Space Telescope Sci¬ 
ence Institute, which is operated by the Association of 
Universities for Research in Astronomy, Inc., for NASA, 
under contract NAS 5-26555. 


APPENDIX 


BUOYANCY RESONANCES 

The inner and outer spiral shocks are due to the steepening of spiral density waves which are excited by the planet 
at Lindblad resonances. At Lindblad resonances, the Doppler shifted frequency matches the disk epicyclic frequency 
{m{Qp — Q) = Azk) and density waves are excited. 

However, besides the epicyclic frequency, the disk also has other natural frequencies. When the disk is not strictly 
isothermal, it has a non-zero Brunt-Vaisala frequency 


N{z) = 


7-1 g{z) 


T ^s,iso 

where = p/p. Matching the Brunt-Vaisala frequency with the Doppler shifted frequency, we have 


7 — I Q.k{R)z 


7 


H 


. X -3/2 

^ ) = ±m{^p - Q). 


(Al) 


(A2) 


Given a m. Equation (A 2 ) gives the position of the resonances. These buoyancy resonances were discovered in 
shearing box simulations (Zhu, Stone, & Rafikov 2012) and studied analytically in Lubow & Zhu (2014). They have 
significant contributions to the planetary torque, especially around the planet, which may affect planet migration. 
These resonances are infinitely thin, and no waves are excited to carry the deposited angular momentum and energy 
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Fig. 21. — Temperature fluctuations (the left panel) and vq (the middle panel) at 2 : = 2H for SMITI. The right panel i s the same as the 
left panel but in Cartesian coordinates. The dotted lines/curves are the position of buoyancy resonances from Equation \A3\ . 
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R 

Fig. 22.— The temperature structure for SM2T 1 at an azimuthal slice (cf) — (/)p = 100®). Spiral shocks and buoyancy resonances are 
labeled. The dotted lines are again from Equation ( |A3| >. 


away. Their dissipation relies on microscopic viscosity or radiative cooling. Thin density ridges with large temperature 
and velocity variations appear at these resonances. 

When various m modes overlap with each other, we can roughly estimate the position of the final density ridges 
caused by buoyancy resonances following Equation 10 and 11 in Zhu et al. (2012). First, given a m, we can calculate 
the corresponding resonance position at R and z. Then using the azimuthal wavelength for this mode A = 27r/m, the 
geometric location of the constant phase 2niT {n is integer) is given by 0 = nA (assuming the phase of buoyancy waves 
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is 0 at the planet position.). Thus, 

(j) = ±2n7r{Qp — ( 2 ) 


H 


7-1 i}K{R)z 


i?2 ) 


(A3) 


We plot tempera ture fluctuations and vq dit z = 2H for SMITI in Fig ure The positions of buoyancy resonances 
given by Equation (A3) are plotted wit h dotted lines/curves. Figure [21] shows that both temperature fluctuations and 
V 0 are nicely tracked by Equation (A3), suggesting that buoyancy resonances exist in disks even with Tcooi = 1- 
Even though buoyancy resonances can affect planet migration, they may not be observed through direct imaging 
technique since the density fluctuations caused by these resonances are much weaker than the spiral density waves 
excited by Lindblad resonances. For example, in the 2 : = 3i7 panels of Figure we can see some density fluctuations 
close to the corotation region (especially for SM2T1, the 1 Mj case), but they are much weaker than the spiral shocks. 

These buoyancy resonances also have vertical structure, as shown in Figure |22]_We sliced through the disk at a fixed 
0, and the buoyancy resonance curves are plotted as dotted curves in Figure 1 2 2| At the disk midplane, = 0 and 
there are no buoyancy resonances. Since N increases with disk height, \ flp — Hjalso needs to increase with height to 
match the Brunt-Vaisala frequency considering A and m are the same at the same 0. Thus, the curves move away 
from the planet position towards the disk atmosphere. 

Figure!^ also shows the inner and outer spiral shocks which are hotter at the shock position. 
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